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ABSTRACT 

Cosmological hydrodynamic simulations of large scale structure in the universe have shown that accretion shocks 
and merger shocks form due to flow motions associated with the gravitational collapse of nonlinear structures. 
Estimated speed and curvature radius of these shocks could be as large as a few 1000 km/s and several Mpc, 
respectively. According to the diffusive shock acceleration theory, populations of cosmic-ray particles can be 
injected and accelerated to very high energy by astrophysical shocks in tenuous plasmas. In order to explore the 
cosmic ray acceleration at the cosmic shocks, we have performed nonlinear numerical simulations of cosmic ray (CR) 
modified shocks with the newly developed CRASH (Cosmic Ray Amr SHock) numerical code. We adopted the 
Bohm diffusion model for CRs, based on the hypothesis that strong Alfven waves are self-generated by streaming 
CRs. The shock formation simulation includes a plasma-physics-based “injection” model that transfers a small 
proportion of the thermal proton flux through the shock into low energy CRs for acceleration there. We found 
that, for strong accretion shocks, CRs can absorb most of shock kinetic energy and the accretion shock speed is 
reduced up to 20 %, compared to pure gas dynamic shocks. For merger shocks with small Mach numbers, however, 
the energy transfer to CRs is only about 10-20 % with an associated CR particle fraction of 10“^. Nonlinear 
feedback due to the CR pressure is insignificant in the latter shocks. Although detailed results depend on models 
for the particle diffusion and injection, these calculations show that cosmic shocks in large scale structure could 
provide acceleration sites of extragalactic cosmic rays of the highest energy. 
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I. INTRODUCTION 

Shocks are ubiquitous in astrophysical environments: 
a few examples are Earth’s bow shock, interplanetary 
shocks, stellar wind terminal shocks, supernova rem¬ 
nants, shocks in radio jets, merger shocks in intracluster 
media, and accretion shocks associated with large scale 
structure formation. Most astrophysical shocks are so- 
called “collisionless shocks” which form in a tenuous 
plasma via electromagnetic “viscosities,” i.e., collective 
electromagnetic interactions between the particles and 
the fields. Hence the magnetic field, especially its irreg¬ 
ular component, is vital to the shock formation process. 
Our discussion will focus on a quasi-parallel shock, in 
which the direction of propagation is almost parallel to 
the magnetic field lines. According to plasma simula¬ 
tions of quasi-parallel shocks (Quest 1988), the particle 
velocity distribution has some residual anisotropy in 
the local fluid frame due to the incomplete isotropiza- 
tion during the collisionless shock formation process 
and so some particles can stream back upstream of 
the shock. Streaming motions of high energy parti¬ 
cles against the background fluid generate strong MHD 
Alfven waves upstream of the shock, which in turn scat¬ 
ter particles and prevent them from escaping upstream 
(e.g., Wentzel 1974; Bell 1978;Quest 1988;Lucek & Bell 
2000). Due to these self-generated MHD waves ther¬ 
mal particles are confined and advected downstream, 
while some suprathermal particles in the high energy 
tail of the Maxwellian velocity distribution may re-cross 


the shock upstream. Then these particles are scattered 
back downstream by those same waves and can be ac¬ 
celerated further to higher energies via Fermi first order 
process. Hence the nonthermal, cosmic-ray particles 
are natural byproducts of the collisionless shock for¬ 
mation process and they are extracted from the shock- 
heated thermal particle distribution (Maikov & Volk 
1998, Maikov & Drury 2001). This “thermal leak¬ 
age” injection process has been observed well in the 
Earth’s bow shock and interplanetary shocks (Ellison, 
Mobius, & Paschmann 1990, Baring et al. 1997). Also 
there have been observational evidence and theoretical 
studies that strongly support the idea that the cosmic 
rays are accelerated via the “diffusive shock acceler¬ 
ation (DSA)” process at various astrophysical shocks 
{e.g., Drury 1983; Blandford & Eichler 1987). 

According to hydrodynamic simulations of large scale 
structure formation {e.g., Kang et al. 1994a; Miniati et 
al. 2000), accretion shocks are formed in the baryonic 
component around non-linear structures collapsed from 
the primordial density inhomogeneities as a result of 
gravitational instability. Those structures can be iden¬ 
tified as pancake-like supergalactic planes, still denser 
filaments, and clusters of galaxies that form at intersec¬ 
tions of pancakes in any variants of the many cosmo¬ 
logical models. These structures are surrounded by the 
hot gas heated by the accretion shocks and the CRs 
(ions and electrons) can be accelerated to very high 
energies at these shocks via the first order Fermi DSA 
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process (Kang, Ryu & Jones 1996, Miniati et al. 2001a, 
b). The accretion shocks around the clusters of galax¬ 
ies could involve flows as fast as a few 1000 kms~^ 
and, so, could be fast acceleration sites for the high en¬ 
ergy cosmic rays up to several x 10^® eV, provided that 
the magnetic field around the clusters is order of mi¬ 
crogauss. Norman, Melrose & Achterberg (1995) also 
suggested that cosmic accretion and merger shocks can 
be good acceleration sites for ultra-hight energy CRs 
above 10^®'^ eV if a primordial held of order ^ 10“®G 
exists, or if microgauss fields can be self-generated in 
shocks. In fact a magnetic field with a typical strength 
of a few microgauss and a principal length scale of 10- 
100 kpc was detected in several clusters of galaxies by 
the Faraday rotation of radiation from distant radio 
sources {e.g., Kim, Kronberg, Tribble 1991, Kronberg 
1994, Taylor et al. 1994, Feretti et al. 1995, Clarke et 
al. 2001, Carilli & Taylor 2002). On the other hand, 
the magnetic fields derived from observed hard X-ray 
and EUV emissions, on the assumption that these are 
inverse Compton (IC) scattering of CMB photons, are 
somewhat lower (~ 0.4 microgauss) (Fusco-Femiano et 
al. 1999, 2000). The origin and evolution of cosmic 
magnetic fields is beyond the scope of this paper. Mag¬ 
netic fields may have been injected into the ICM by ra¬ 
dio galaxies. They may have been seeded at shocks in 
the course of structure formation, and then stretched 
and amplified up to 0.1-1 microgauss levels by turbulent 
flow motions in ICM and also in filaments and sheets 
(Kulsrud et al. 1997, Ryu, Kang & Biermann 1998). 

Recent observations in EUV and hard X-ray have re¬ 
vealed that some clusters possess excess radiation com¬ 
pared to what is expected from the hot, thermal X- 
ray emitting ICM {e.g., Sarazin & Lieu 1998; Lieu et 
al. 1999; Ensslin et al. 1999; Fusco-Femiano et al. 1999; 
Sarazin 1999). One mechanism proposed for the ori¬ 
gin of this component is the IC scattering of cosmic 
microwave background photons by CR electrons accel¬ 
erated by merger shocks and accretion shocks around 
the clusters. Also it has been suggested that the diffuse 
gamma-ray background radiation could originate from 
the same process (Loeb & Waxman 2000, Miniati 2002, 
Scharf & Mukherjee 2002). The same mechanisms that 
are capable of producing CR electrons may have pro¬ 
duced CR protons, although the existence of CR pro¬ 
tons in the ICM has not yet been directly observed. 
The existing evidence for substantial CR populations 
in these environments argues that nonthermal activi¬ 
ties in the ICM could be important in understanding 
the dynamical status and the evolution of clusters of 
galaxies (Sarazin & Lieu 1998; Lieu et al. 1999). CR 
protons and electrons may provide a significant pres¬ 
sure to the ICM, perhaps, comparable to the thermal 
gas pressure (Lieu et al. 1999, Colafrancesco 1999), as 
it is for the galactic CRs in the ISM of our galaxy. Col¬ 
lisions of CR protons in the ICM generate a flux of 
7 -ray photons through the production and subsequent 
decay of neutral pions. While such y-rays have not 
yet been detected from clusters recent estimates have 


shown that y-ray fluxes from the nearest rich clusters, 
such as Coma, are within the range of what may be 
detected by the next generation of y-ray observatories 
(Ensslin et al. 1997, Sreekumar et al. 1996, Miniati et 
al. 2001). 

According to DSA theory a significant fraction (up 
to 90%) of the kinetic energy of the bulk flow associ¬ 
ated with the strong shock can be converted into CR 
protons, depending the CR injection rate (Drury 1983, 
Jones & Kang 1990, Berezhko, Ksenofontov, & Yelshi 
1995). If as much as lO”"* — 10“® of the particle flux 
passing through the shock were injected into the CR 
population, the CR pressure would dominate and the 
nonlinear feedback to the underlying flow would be¬ 
come substantial. Recently Gieseler, Jones & Kang 
(2000) have developed a novel numerical scheme that 
self-consistently incorporates the thermal leakage in¬ 
jection based on the analytic, nonlinear calculations 
of Maikov (1998). This injection scheme, which has 
only one tightly restricted adjustable parameters, has 
been implemented into the combined gas dynamics and 
the CR diffusion-convection code with the Adaptive 
Mesh Refinement technique by Kang, Jones & Giesler 
(2002). The CR injection and acceleration efficien¬ 
cies at quasi-parallel, plane-parallel shocks were calcu¬ 
lated with their new numerical code named as CRASH 
(Comic-Ray Amr SHock) code. They found that about 
10“® of incoming thermal particles are injected into the 
CRs, that up to 60 % of initial shock kinetic energy 
is transferred to CRs for strong shocks, and that the 
shock speed is reduced up to ~ 17 % for shocks with 
Mach number greater 30. These results have confirmed 
the findings of previous studies which adopted a simpler 
injection model that included a fully adjustable free 
parameter {e.g., Berezhko et al. 1995, Kang & Jones 
1995). 

In this contribution we will present the numerical 
simulation results for quasi-parallel shocks in ID plane- 
parallel geometry with the physical parameters rele¬ 
vant for the cosmological shocks emerging in large scale 
structure formation of the Universe. In the next sec¬ 
tion the details of numerical simulations, including the 
basic equations, numerical method, and model parame¬ 
ters, will be given. The simulation results are presented 
and discussed in §III, followed by a brief summary in 
§IV. 

II. Numerical Methods 
(a) Basic Equations 

We solve the standard gasdynamic equations with 
GR pressure terms added in the conservative, Eulerian 
formulation for one dimensional plane-parallel geome¬ 
try: 


^ d{up) 
dt dx 


(1) 
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d{pu) d{pv? + Pg + Pc) 
dt dx 


( 2 ) 


d{peg) djpegU + PgU + Pcu) 

dt dx 


-L{x,t), 


( 3 ) 


^ djuS) 
dt dx 


( 4 ) 


where Pg and are the gas and the CR pressure, re¬ 
spectively, Cg = Pg/lpijg — 1)] -I- u^/2 is the total en¬ 
ergy density of the gas per unit mass and the rest of 
the variables have their usual meanings. The injection 
energy loss term, L{x, t), accounts for the energy of the 
suprathermal particles injected to the CR component 
at the subshock. Here, S = Pg/p'^‘>~^, the “Modified 
Entropy”, is introduced in order to follow the preshock 
adiabatic compression accurately in strong CR mod¬ 
ified shocks. We note that the equation (4) is valid 
only outside the dissipative subshock; i.e., where the 
gas entropy is conserved. Hence the modified entropy 
equation is solved outside the subshock, while the total 
energy equation is applied across the subshock. 

The diffusion-convection equation for the pitch angle 
averaged CR distribution function, f{p, x, t), e.g., Skilling 
1975) is given by 


df df 1 du df d df 


and k{x,p) is the spatial diffusion coefficient. For con¬ 
venience we always express the particle momentum, p 
in units rUpC. As in our previous studies, the function 
g{p) = p'^f{p) is evolved instead of /(p) and y = ln{p) 
is used instead of the momentum variable, p for that 
step. Then the CR pressure is calculated from the non- 
thermal particle distribution as follows: 

. 4 2 r . 3 dp 

Pc =-TrmpC / g{p)—===. (6) 

a Jo VP + 1 


(b) CRASH: CR/AMR Hydrodynamics Code 


diffusion transport model with a steeply momentum- 
dependent diffusion coefficient, the highest energy, rel¬ 
ativistic particles have diffusion lengths many orders 
of magnitude greater than those of the lowest energy 
particles. 

To follow the acceleration of highly relativistic CRs 
from suprathermal energies, all those scales need to be 
resolved numerically. However, the diffusion and accel¬ 
eration of the low energy particles are important only 
close to the shock owing to their small diffusion lengths. 

Thus it is necessary to resolve numerically the diffusion 
length of the particles only around the shock. To solve 
this problem generally we have developed the CRASH 
code by combining a powerful “Adaptive Mesh Refine¬ 
ment” (AMR) technique (Berger & Le Veque 1998) and 
a “shock tracking” technique (Le Veque & Shyue 1995), 
and implemented them into a hydro/CR code based 
on the wave-propagation method (Kang et al. 2001; 
Kang et al. 2002). The AMR technique allows us to 
“zoom in” inside the precursor structure with a hier¬ 
archy of small, refined grid levels applied around the 
shock. The shock tracking technique follows hydro- 
dynamical shocks within regular zones and maintains 
them as true discontinuities, thus allowing us to refine 
the region around the gas subshock at an arbitrarily 
fine level. The result is an enormous savings in both 
computational time and data storage over what would 
be required to solve the problem using more traditional 
methods on a single fine grid. 

(c) Injection Model 

In the “thermal leakage” injection model, some suprather¬ 
mal particles in the tail of the Maxwellian distribution 
swim successfully against the Alfven waves advecting 
downstream, and then leak upstream across the sub¬ 
shock and get injected in the CR population. In order 
to model this injection process in Gieseler et al. (2001) 
we adopted a “transparency function”, Tesc, which ex¬ 
presses the probability that supra-thermal particles at 
a given velocity can leak upstream through the mag¬ 
netic waves, based on non-linear particle interactions 
with self- generated waves (Maikov and Volk 1998). In 
this scheme, the transparency function is approximated 
by the following functional form, 


Unlike ordinary gas shocks, the CR shock includes 
a wide range of length scales associated not only with 
the dissipation into “thermal plasma”, but also with 
the nonthermal particle diffusion process. Those are 
characterized by the so-called diffusion lengths, 

L'diff(p) = n(jp)/u, (7) 

where k{p) is the spatial diffusion coefhcient for CRs of 
momentum p, and u is the characteristic flow velocity 
(Kang & Jones 1991). Accurate solutions to the CR 
diffusion-convection equation require a computational 
grid spacing significantly smaller than Udiff, typically. 
Ax ~ O.llldiff(p)- On the other hand, for a realistic 


Tesc{e,v/ud) = H[v - {1 + e)]{l - ^(1-4) 

•exp{-[b-(l + e)]-2}, (8) 

which depends on the ratio of particle velocity, v, to 
downstream flow velocity in the subshock rest-frame, 
Ud. Here PI is the Heaviside step function. The in¬ 
verse wave-amplitude parameter, e = is de¬ 

fined in Maikov and Volk (1998) and measures the ra¬ 
tio of the amplitude of the postshock MHD wave tur¬ 
bulence B^ to the general magnetic field aligned with 
the shock normal, Bq. Here v = e v/ud is the nor¬ 
malized particle velocity. This function behaves like a 
smoothed step function, and the injection takes place in 
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momentum space where the dTgsc/dp is greatest. The 
breadth of the thermal velocity distribution relative the 
downstream flow velocity in the subshock rest-frame 
(i.e., Uth/wd) determines the probability of leakage, 
and so the injection process is sensitive to the veloc¬ 
ity jump at the subshock, which depends on the sub¬ 
shock Mach number. The injection rate increases with 
the subshock Mach number, but becomes independent 
of Ms in the strong shock limit of Mg > 10 (Kang et 
al. 2002). Thus the injection rate should change with 
the shock strength in the time-dependent evolution of 
CR modified shocks, as the subshock weakens by way 
of precursor compression. 

The only free parameter of the adopted transparency 
function is e and it is rather well constrained, since 
0.3 ^ e ^ 0.4 is indicated for strong shocks (Maikov 
& Volk 1998). However, it turns out the injection rate 
depends sensitively on the value of e, due to the ex¬ 
ponential cut off in a thermal velocity distribution. It 
is also expected that the wave generation is weaker for 
low Mach shocks, leading to the larger values of e. So 
in this study we will consider 0.2 < e < 0.4. 

In the CRASH code, in order to emulate numerically 
thermal leakage injection, we first estimate the number 
of suprathermal particles that cross the shock according 
to the diffusion-convection equation, and then we allow 
only a small fraction of the combined advective and 
diffusive fluxes to leak upstream with the probability 
prescribed by Tesc- The readers are referred to Kang et 
al. (2002) for more details of our numerical scheme for 
thermal leakage injection model. 

III. Numerical Model 

(a) Diffusion Coefficient Model 

Diffusive acceleration at shocks depends on the ex¬ 
istence of Alfvenic turbulence capable of strongly scat¬ 
tering energetic protons. The irregularities in the mag¬ 
netic field within the cluster can be generated by var¬ 
ious dynamical effects such as the supersonic motion 
of galaxies through ICM, mergers of substructures and 
galactic winds. In addition, according to the cosmologi¬ 
cal hydrodynamic simulations, the downstream regions 
of the accretion shocks (be., ICM) are full of turbulent 
flow motions. Some observations (Feretti et al. 1995) 
indicate that the ICM field is likely to be tangled on 
scales of the order of less than 1 kpc. Then turbulent 
flows may subsequently generate the irregularities in 
magnetic field. Outside the accretion shock, however, 
we must hypothesize the existence of pre-existing field 
turbulence. Lucek & Bell (2000) have shown by nu¬ 
merical simulations that CR streaming induces large- 
amplitude Alfven waves in a quasi-parallel shock, im¬ 
plying that the shock formation process itself can gen¬ 
erate the necessary turbulent field. So it is commonly 
assumed in calculations of diffusive shock acceleration 
that the scattering Alfven waves are self-generated by 
the CRs themselves by the so-called “streaming insta¬ 


bility” (Wentzel 1974). Adiabatic compression of that 
field and turbulence by convergence of the inflow onto 
the accretion shock might very reasonably lead to an 
acceptable level of magnetic irregularities in the shock 
vicinity. 

The Bohm diffusion model represents a saturated 
wave spectrum and gives the minimum diffusion coeffi¬ 
cient as Kb = l/3rgU, when the particles scatter within 
one gyration radius (rg) due to random scatterings off 
the self-generated waves. This gives 

AC(P) =Ko ^^2^1)1/2 ’ (9) 

where Ko = 3.13 x 10^^cm^s“^i?“^ and is the mag¬ 
netic field strength in units of microgauss. In order to 
model amplification of self-generated turbulent waves 
due to compression of the perpendicular component of 
the magnetic field, the spatial dependence of the diffu¬ 
sion is modeled as 

k(x,p) = k(p)(-^), (10) 

p{x) 

where po is the upstream gas density. This form is 
also required to prevent the acoustic instability of the 
precursor (Kang, Jones & Ryu 1992). 


(b) Time and Length Scales 


The mean acceleration time scale for a particle to 
reach momentum p is determined by the velocity jump 
at the shock and the diffusion coefhcient (e.g., Drury 
1983), that is, 


<dp > 
^ dt ^ 


-(^ + ^). (11) 

Ml — M2 Ml M2 


Here the subscripts, 1 and 2, designate the upstream 
and downstream conditions, respectively. If the strong 
shock limit is taken (be., U 2 = mi/ 4), and we assume for 
a turbulent field that B/p is constant across the shock 
(be., k/u = constant), then Tacc ~ Skb/u^- Using the 
Bohm diffusion coefficient in equation (9), the mean 
acceleration time scale is given by 


( 12 ) 

This shows that the diffusive acceferation time scale 
increases directly with the particle energy and that it 
takes about a billion years to accelerate the proton to 
E = 10^®eV at a typical cluster shock, assuming ~ 
1. The diffusion length of the CR protons is given by 


Ddiff 


Kb 

Us 


1.02Mpc( 


IQio 




1000km 


-1 


(13) 

So the CR protons oi E = 10^®eV diffuse on the length 
scale comparable to the cluster size. Note the mildly 
relativistic protons (p ~ 1 — 10) are almost instanta¬ 
neously accelerated at these shocks compared to the 
cosmological time scale and diffuse on the length scale 
much smaller than the cosmological scale. 
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Fig. 1 .— Time evolution of the shock driven by an accretion flow with po = 1, wo = — 1 and Mq — 30 which is reflected at 
a; = 0. Seven levels of refinements (Imax = 7) were used and the inverse wave-amplitude parameter e = 0.2 was adopted. The 
snapshots are shown at t = 5, 10, 50, 150, and 200. The shock propagates to the right, so the leftmost profile corresponds 
to the earliest time. For t — 200, data at each cell is shown as filled circles to show clearly the subshock jump. Note the 
distance from the reflecting plane is in a logarithmic scale. 


(c) ID Plane-parallel Shock Models 

In this contribution, we study the CR accelera¬ 
tion and its dynamical effects at one-dimensional (ID) 
quasi-parallel shocks which form due to the accretion 
flows in a plane-parallel geometry. In general, cosmic 
shocks associated with large scale structure formation 
can be oblique and have various geometries, depending 
on the types of nonlinear structure onto which the ac¬ 
cretion flow falls. Roughly speaking we have: ID plane- 
parallel shocks around the sheets, 2D cylindrical shocks 
around the filaments, and 3D spherical shocks around 
the clusters of galaxies. Due to severe requirements on 
the computational resources, our simulations can follow 
the acceleration of the protons from suprathermal ener¬ 
gies (p ~ 10“^) to mildly relativistic energies (p ~ 50). 
For the particles in this energy range, the acceleration 
time scales (< 40 years) are much shorter than and the 


cosmological time scale (tn 10^° years) and the diffu¬ 
sion length scales (< O.lpc) are much smaller than the 
curvature of multi-dimensional cosmic shocks {Rs ~ a 
few Mpc). On those scales where D^is ^ Rs, the dif¬ 
fusion and acceleration of CRs can be studied with the 
ID plane-parallel shock models. Thus, for our purposes 
the shocks that occur in major mergers of substructures 
can be represented approximately by ID plane-parallel 
shocks. As shown in Kang et al. (2002), the maxi¬ 
mum CR pressure and the shock structure modifica¬ 
tion approach to time asymptotic limits as the parti¬ 
cles become relativistic p ~ 1. This is due to the bal¬ 
ance between acceleration and diffusion of CRs, that is, 
fewer particles are accelerated to higher energies, and 
higher energy CRs diffuse over larger volume of space 
as the CR shock evolves. Thus our simulations can pre¬ 
dict long-term behaviors of CR modified shocks, even 
though our integration time is only a small fraction of 
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Fig. 2 .— Same as Figure 1 except Mq = 5. 


tn- 

Here we ignore gravity, so changes in the infall veloc¬ 
ity and density, and the background cosmological evolu¬ 
tion, because the integration time of our simulations is 
much shorter than any of the cosmological time scales. 
So we consider a pancake shock formed by the steady 
accretion flow with a constant density and pressure: a 
ID simulation box with ^^Xmax] and an accretion flow 
entering into the right boundary of the simulation box 
with a constant density, pq, pressure, Pg, 0 ) and velocity, 
Mq. The flow is reflected at the left boundary (x = 0) 
(i.e., pancake middle plane). A shock forms and prop¬ 
agates to the right. For a hydrodynamic shock with¬ 
out the CRs, the shock speed is = |Mo|/(r — 1) in 
the simulation frame and Ug = |uo|r/(r — 1) in the far 
upstream rest frame, where r is the compression ra¬ 
tio across the shock. Throughout the paper we denote 
the velocities in the simulation frame by the unprimed 
variables, while the velocities in the rest frame of far 
upstream flow by the primed variables. For strong hy¬ 
drodynamics shocks, r = 4 and u'g = (4/3)|mo|- A CR 


modified shock consists of a smooth precursor and a 
subshock, since the CRs diffuse upstream of the sub¬ 
shock, and the CR pressure decelerates and heats the 
preshock flow adiabatically, resulting in weakening of 
the subshock. We denote the values immediately up¬ 
stream of the subshock as ui, pi, and and the 
values immediately downstream of the subshock as U 2 , 
P 2 and Pg^ 2 - We also denote the subshock speed rela¬ 
tive to the immediate upstream flow as rtsub, which is 
smaller than u'g due to the precursor deceleration. 

We designate the shock models by the Mach number 
of the accretion flow that is defined by 


Mq = 


kol 

(7-fg.o/po)^A’ 


(14) 


where once again uq, po, Pg^ are the accretion velocity, 
the density and pressure of the infall flow, respectively. 
The Mach number of the accretion shock resulting from 
such accretion flow is Mg « 4/3Mo for strong gasdy- 
namic shocks. We set the far upstream density and flow 
values as po = 1) Wo = — 1 in code units for all models. 
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11^=2. £=0.S, 1„=7, t= 5. 10, 50, 
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Fig. 3 .— Same as Figure 1 except Mo = 2. 


while we vary the gas pressure, Pg,o for different values 
of Mq, where 2 < Mq <100. So the preshock gas is 
colder for high values of Mq. 

The following parameters are used: the gas adiabatic 
index, jg = 5/3 and the velocity normalization con¬ 
stant Uo = 1500 kms“^ {(3 = Uo/c = 0.005). Then the 
diffusion coefficient, Kq = 3.13 x 10^^cm^s“^R“^, de¬ 
termines the length and time scales as Xo = Ko/uo and 
to = Ko/Uqi respectively, which represent, in fact, the 
diffusion length and time scales for the protons of p = 1. 
For Xo = 2.1 X lO^"* cm and to = 1.4 x 10®s. 

These scales are much smaller than cosmological length 
and time scales as discussed before, which justifies 
our assumptions about the steady accretion flow and 
the ID plane-parallel geometry. The gas density and 
pressure normalization constants, po and Po = PoU^, 
are arbitrary. For the hydrogen number density of 
nn = 10“^cm“^, for example, po = 2.34x 10“^®g cm“^ 
and Po = 5.26 x 10“^^erg cm“^. So if one assumes the 
Bohm diffusion model, the three parameters, i.e., (3, 
Rp, and Mq, sufficiently define a CR shock model. 


100. 150. 200 



log(xj 



iog(x) 


Throughout the paper and in the code physical vari¬ 
ables are given in units of the normalization constants, 
Xo, to, Uo, Po, and Po- 

Although the theoretically preferred values of the 
inverse wave-amplitude parameter, e, lie between 0.3 
and 0.4 for strong shocks (Maikov 1998), such values 
lead to very efficient initial injection and most of the 
shock energy is transfered to the CR component for 
strong shocks of Mq ^ 30 (Kang et al. 2002). As a 
more conservative option we have considered a set of 
models for 2 < Mq < 100 with e = 0.2 and another 
set of models for Mq = 2 with 0.2 < e < 0.4. The 
former is chosen to explore the dependence of the CR 
acceleration on the accretion flow Mach number for a 
given value of e, while the latter is chosen to explore 
the dependence on the injection rate for a low Mach 
number flow. 

The simulations were carried out on a base grid with 
Axq = 0.02 using = 7 additional grid levels, so 
Axr = 1.56 X 10“^ on the finest grid. The simulated 
space is a; = [0, 200] and N = 10000 zones are used 
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M„=2, €=0.4, 1„=7, t= 5. 10, 50, 50, 150, 800 



log(xj log(xj 
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Fig. 4. — Same as Figure 1 except Mq = 2 and e = 0.4. 


on the base grid for the models that were integrated 
to t/to = 200. For the models that were integrated to 
t/to ^ 100, X = [0,100] and N = 5000 are used. The 
number of refined zones around the shock is N^f = 100 
on the base grid and so there are 2Nrf = 200 zones on 
each refined level. The length of the refined region at 
the base grid is 2, so 1/100 of the entire simulated space 
on the base grid is refined. To avoid technical difficul¬ 
ties, the multi-level grids are used only after the shock 
propagates away from the left boundary at the distance 
of = 1. With < 1, the downstream refined re¬ 
gion is outside the left boundary of the simulation box 
and the full length of the refined region around the 
shock cannot be set down with the current version of 
the CRASH code. After the shock moves to a;* = 1 (at 
t « 3 for strong shocks), the AMR technique is used 
and the CR injection and acceleration are activated. 
This initial delay of the CR injection and acceleration 
should not affect the final outcomes. We integrate the 
simulation until t = 100 — 200, so that the maximum 
momentum achieved by the end of simulation is of order 


of Pmax ~ 40, above which the CR distribution func¬ 
tion decreases exponentially. For all models we use 230 
uniformly spaced logarithmic momentum zones in the 
interval log(p/mpc) = [logpo, logpi] = [-3.0,-1-3.0] 

For gasdynamic variables the continuous boundary 
condition is used at right boundary, while the reflecting 
boundary condition is applied to left boundary of the 
simulation box. For the CR distribution function a 
continuous boundary is assumed for the advection step 
and a no-flux boundary condition is adopted for the 
spatial diffusion step. Either below po or above pi; 
g{p) = 0 is assumed. 

IV. RESULTS 

(a) Modified Shock Structure 

We show the time evolution of shocks at t = 5, 
10, 50, 100, 150 and 200 for models with the accre¬ 
tion flow Mach number, Mq = 30, 5, and 2 in Fig¬ 
ures 1-3. The inverse wave amplitude parameter was 
assumed to be e = 0.2. As CRs are injected and ac- 
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celerated at the shock, the CR pressure increases and 
diffuses upstream, leading to a precursor in which the 
upstream flow is decelerated and compressed adiabati- 
cally. As the CR precursor grows, the subshock slows 
down and the postshock density increases, while the 
postshock gas pressure decreases. Because the injection 
rate is quite high for strong shocks, the CR energy in¬ 
creases and the modification to the flow structure pro¬ 
ceeds very quickly in a time scale comparable to face 
for p ~ 2 — 3. 

Kang et al. (2002) also considered similar shock 
models, but with different flow conditions at the left 
boundary. There a gasdynamic shock was set up ini¬ 
tially and allowed to evolve with open boundaries both 
far upstream and downstream. Because the CRs are 
not allowed to diffuse out at the reflecting, downstream 
boundary in the current models, the CR pressure builds 
up faster. Otherwise, however, the CR shocks evolve 
similarly in both models: 1) The total transition con¬ 
sists of a precursor and a subshock that weakens to 
a lower Mach number shock, but does not disappear 
entirely. 2) After an initial quick adjustment, the 
CR pressure at the shock reaches approximate time- 
asymptotic values when the fresh injection and accel¬ 
eration are balanced with advection and spreading of 
high energy particles due to strong diffusion. For strong 
shocks, T’c,2/po(Ms o)^ ^ 0.56, where u'^ q is the shock 
speed before any significant nonlinear CR feedback oc¬ 
curs. 3) Once the postshock CR pressure becomes con¬ 
stant, the shock structure evolves approximately in a 
“self-similar” way, because the scale length of shock 
broadening increases linearly with time. 4) A postshock 
“density spike” forms due to the nonlinear feedback of 
the CR pressure (see also Jones & Kang 1990). 5) For a 
given inverse wave-amplitude parameter, e, the CR ac¬ 
celeration efficiency and the flow modification depend 
sensitively on the shock Mach number, but they seem 
to converge at the strong shock limit {Mg > 30). 

For the Mq = 30 model, for example, the initial 
unmodified shock speed is q = 4/3 (Mach num¬ 
ber Mg = 40), but the shock slows down to u'g « 
1.1 due to CR nonlinear modification. So the ratio 
Pc, 2 /Poiu'g)^ 0.8, where u'g is the instantaneous 

shock speed. For strong shocks the CR pressure dom¬ 
inates over the gas pressure with a CR injection frac¬ 
tion ~ 10“^ — 10“^ (see Figure 8 below), and a ra¬ 
tio Fc. 2 /po(m's)^ ^ 0.8 — 1.0, which is consistent with 
the simulation results of Kang et al. (2002). For the 
Mq = 5 model, Pc, 2 /Poiu'g 0.27 and still shows 

substantial degree of the flow modification. For Mq = 2 
model, however, the CR pressures increases only to 
Pc,2/Poiu'sY ~ 0-03 and the flow structure is not mod¬ 
ified at all. According to Miniati et al. (2000), the cos¬ 
mic shocks inside the ICM {i.e., merger shocks) have 
mostly low Mach numbers of 1 ^ Mg ^ 5, because 
these shocks propagate into the already hot medium. 
But the shock speed is not much different from the ac¬ 
cretion shocks that propagate into the cold medium and 
so have high Mach numbers. Thus the sample models 


shown in Figures 1-3 can provide qualitative pictures 
on the effects of the CR acceleration to dynamics of 
cosmological shocks. 

We note, however, that the value of e = 0.2 is per¬ 
haps unrealistically small for a low Mach number shock, 
since the self-generation of waves is expected to be less 
efficient. So, we considered three additional models for 
Mq = 2 (Mg = 3) with larger values of e = 0.25, 0.3, 
and 0.4. The CR pressure increases with higher values 
of e due to higher leakage fractions, as expected. So, 
for example, the ratio Pc. 2 /po('w/)^ ~ 0.22 for e = 0.4 
model (see Figure 4). But even in this rather high in¬ 
jection model with a CR fraction of ^ ~ 10“^, the flow 
modification is minimal except for the reduction of the 
gas pressure. As expected, the CR acceleration is very 
inefficient and the CR nonlinear feedback is insignifi¬ 
cant in weak shocks of a Mach number of a few. 

(b) Particle Distribution Function 

Evolution of the CR distribution function at the 
shock is given for the models of Mq = 30,10, 5, and 
2 with e = 0.2 in Figure 5. Thermal particles are rep¬ 
resented by a Maxwellian distribution below p ~ 10“^. 
One can see that the Maxwellian distribution shifts to 
lower momenta, as the postshock temperature reduces 
due to energy transfer to CRs. Just above the injection 
pool, the distribution function changes smoothly from 
the thermal distribution to an approximate power-law 
whose index is close to the test-particle slope for the 
subshock. While a fraction of particles injected ear¬ 
lier continues to be accelerated to higher momenta, so 
that Pmax(0 increases, the amplitude of g{p) at the 
shock and at a given momentum decreases with time 
for p < Pm&xit) due to diffusion. The distribution func¬ 
tion g{p) shows the characteristic “concave upwards” 
curves reflecting modified shock structure (including 
the precursor) for the shocks with Mq > 10. For the 
models with Mq = 2 and e = 0.2, the CR modifica¬ 
tion to the flow structure is insignificant, so the particle 
distribution is a power-law with the test-particle slope 
{i.e., f{p) = /oP~* and q = 3r/(r — 1) w 4.5, where 
r = 3 for Mg = 3). 

In Figure 6, the CR distribution function at the 
shock is shown for the Mq = 2 model for e = 0.2, 
0.25, 0.3, and 0.4. For weaker wave fields (larger e) the 
particles with lower momenta can leak upstream, so the 
injection pool is located at the lower momenta closer 
to the peak of the Maxwellian distribution, resulting 
in a higher injection rate. Even though about 20 % of 
the shock kinetic energy is transfered to the CR com¬ 
ponent for e = 0.4, the flow velocity structure is only 
slightly affected. So the CR distribution function is 
dictated by the test-particle like power-law spectrum 
for Mq = 2 models, roughly independent of the value 
of e. 
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log(p/mpC) log(p/mpC) 

Fig. 5 .— Evolution of the CR distribution function at the shock, represented as <; = is plotted for the models of 

Mo = 30, 10, 5, and 2 with e = 0.2. 


(c) Subshock Evolution 

The speed of the “initial hydrodynamic shock”, u'^ g, 
that emerges from the reflecting plane before the CR in¬ 
jection/acceleration begins depends on the Mach num¬ 
ber of the accretion flow. It is given by g = r/(r — 
l)|uo|, so that u'g Q = 1.5 for Mq = 2, u'^ q = 1.36 for 
Mg = 5, and u'^ q = 4/3 for Mg > 10. After the CR 
injection/acceleration starts, the CR pressure modifies 
the shock flow and the “instantaneous shock” speed 
relative to the far upstream flow, u'^, decreases. We 
plotted the position of the shock, Xs, against time for 
models with e = 0.2 in the upper left panel of Figure 7. 
We see that the shock is decelerated further in higher 
Mach number models, since the CR pressure is greater 
and so is its dynamical feedback to the flow. 

The injection rate in the thermal leakage injection 
model depends on the strength of the subshock, so the 
subshock speed relative to the immediate preshock gas 
in the precursor, Msub, is important. Because of the 


pre-deceleration in the precursor, especially for strong 
shocks, Usub can be much smaller than the shock speed 
relative to the far upstream. We plotted this subshock 
speed in the upper right panel of Figure 7. For Mg = 2 
Msub ~ u'g, but Msub 'C m/ for Mg > 30. We also plotted 
the pre-subshock gas density, pijpo, and the immedi¬ 
ate postshock gas density, pijpo-, in the lower two pan¬ 
els. The modification to the flow velocity occurs mostly 
before t ^ 10, so the subshock velocity decreases at 
the same time, but the precursor compression contin¬ 
ues even after t > 10, especially in the strong shocks. 
We note that the subshock persists and so a completely 
smooth transition never develops by the termination 
time of our simulations. 

(d) Injection and Acceleration Efficiencies 

We define the injection efficiency as the fraction of 
particles that have entered the shock from far upstream 
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Fig. 6 .— Evolution of the CR distribution function at the shock, represented as g = is plotted for the Mq = 2 

models with e = 0.2, 0.25, 0.3 and 0.4. 


and then injected into the CR distribution: 


box from far upstream. 


m 


^^JPl t)p^dp 

ll nou'^{t')dt' 


(15) 




dxEcR(a^,0 

0.5po«,o)^^ 


(16) 


where /cr is the CR distribution function, no is the 
particle number density far upstream, and ti is the time 
when the CR injection/acceleration is turned on (ti « 3 
for Mo > 5 and t « 2 for Mo = 2). If the subshock 
becomes steady, then nou'g is the same as the particle 
flux swept by the subshock, niUsub, where ni is the 
particle number density immediately upstream to the 
subshock. In our simulations these two fluxes can differ 
up to 10 % for strong shock models, because the shock 
speed is changing slowly. 

As a measure of acceleration efficiency, we define the 
“CR energy ratio”; namely the ratio of the total CR en¬ 
ergy within the simulation box to the kinetic energy in 
the initial shock frame that has entered the simulation 


Since our shock models have the same accretion density 
and velocity, but different gas pressure depending on 
Mo, we use the kinetic energy flux rather than the total 
energy flux to normalize the “CR energy ratio”. 

In Figure 8 we show the CR energy ratio, $, the 
CR pressure at the shock normalized to the ramp pres¬ 
sure of the upstream flow in the instantaneous shock 
frame, Pc, 2 /Po{u's)^, and the “time-averaged” injection 
efficiencies, for shocks with different Mach numbers 
when e = 0.2 (left three panels). For all Mach num¬ 
bers the postshock Pc ,2 increases until a balance be¬ 
tween injection/acceleration and advection/diffusion of 
CRs is achieved, and then stays at a steady value af¬ 
terwards. The time-asymptotic value of the CR pres¬ 
sure becomes, once again, Pc.2/po(Rs,o)^ ~ 0-56 and 
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Fig. 7 .— Shock location from the reflecting plane, Xs, the shock velocity relative to the immediate preshock flow, Msub, 
preshock density, pijpo, and postshock density, P 2 /P 0 are plotted as a function of time for Mo = 2 — 100 and t = 0.2. 


Pc ,2 / poiu's)'^ ~ 0.8, for Mq = 30 with e = 0.2. We note 
that the “undulating” features in the time evolution of 
Pc ,2 seem to be numerical artifacts and not real. Un¬ 
like other spatially averaged or integrated quantities, 
the plotted Pc ,2 is sampled exactly at the subshock 
{i.e., from one zone) whose properties show small, noisy 
variations in time. They seem in particular to corre¬ 
spond to times when the subshock crosses a regular 
zone boundary and are most prominent for high Mach 
number models of Mq = 30, and 100. These features 
can be also seen in the preshock and postshock densities 
shown in Figure 7. 

The CR energy ratio, <&, increases with time as 
CRs are injected and accelerated, but it asymptotes 
to a constant value, once Pc ,2 has reached a quasi¬ 
steady value. This results from the approximate “self¬ 
similar” evolution of the Pc spatial distribution. Time- 


asymptotic values of increase with Mq and w 0.5 
for Mq = 30 at the terminal time. 

As discussed in §11. (c) in the thermal leakage model 
the injection rate is higher for higher subshock Mach 
number, because the ratio of the thermal peak velocity 
to the injection velocity is smaller. This ratio, however, 
becomes independent of Mgub for M^ub > 10, so the 
injection rate also shows the same trend. This Mach 
number dependence of ^ can be observed in Figure 8, 
where the initial value of ^ (at t = 3) increases with Mq, 
but it is about the same for Mq > 10. After the ini¬ 
tial quick increase, it decreases in time as the subshock 
weakens due to the pre-deceleration in the precursor. 

Once again, in order to explore the dependence of 
our injection model on the parameter e, especially for 
low Mach shocks, we also show the results for the 
Mq = 2 {Ms = 3) models with e = 0.2, 0.25, 0.3, 
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Fig. 8. — The ratio of total CR energy in the simulation box to the kinetic energy in the initial shock rest frame that has 
entered the simulation box from upstream, $(t), the postshock CR pressure in units of far upstream ram pressure in the 
instantaneous shock frame, and time-averaged injection efficiency, ^(t). Left three panels are for Mo — 2 — 100 and e = 0.2. 
Right three panels show the same quantities for Mo = 2 and e = 0.2, 0.25, 0.3, and 0.4. 


and 0.4 in the right three panels of Figure 8. As ex¬ 
pected, the injection rate is higher for larger values of 
e, so the CR pressure and $ are higher. In Kang et 
al. (2002), we made a similar comparison for a wider 
range of Mach numbers and found that the dependence 
on e is much weaker for stronger shocks. In a physical 
model the parameter, e, should depend on the subshock 
Mach number. But it is not well understood how e 
should vary with the subshock Mach number for weak 
shocks. For strong shocks, the average injection rate is 
about 10“^ with e = 0.2, which corresponds to strong 
wave generation and inefficient leakage. This injection 
rate is in fact in a good agreement with what has been 
observed in the Earth’s bow shock (Quest 1988). For 
Ms = 3 shock, the similar injection rate is obtained for 
0.25 < e < 0.3, and the CR pressure is about 10-15 % 


of the shock ram pressure, which could be considered 
substantial. Then we conclude that CRs can absorb a 
significant portion of the shock kinetic energy at cos¬ 
mological shocks, if about ~ 10“^ of the particles are 
injected into the CR component regardless of the de¬ 
tails of the injection process. 

V. SUMMARY 

We have calculated CR acceleration at ID quasi¬ 
parallel shocks by using our CRASH (Cosmic-Ray Amr 
SHock) code (Kang et al. 2002), which incorporates the 
“thermal leakage” injection process to the CR/hydro 
code that solves the CR diffusion-convection equation 
along with CR modified gasdynamic equations. Our 
simulations are performed in a ID plane-parallel space 
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in which shocks are generated by the accretion flow 
reflecting of! the central symmetry plane. For conve¬ 
nience, we considered the accretion velocity of face = 
1500kms“^ and the magnetic held of 1 microgauss as 
fiducial values for our simulations. However, the gen¬ 
eral conclusions can be applied to similar shock speeds, 
because the CR acceleration is controlled mainly by two 
physical parameters, the shock Mach number, Mg, and 
the inverse wave-amplitude parameter, e. The current 
simulations are similar to those presented in Kang et 
al. (2002), in which freely propagating shocks with open 
boundaries were considered, except that we adopted 
a reflecting boundary in the downstream region and 
our shock velocities, Vg = 2000 — 2250kms“^, are lower 
than their value of 3000 kms“^. Having the reflecting 
plane expedites the CR acceleration in our simulations, 
because CRs are trapped between the shock and the 
downstream, reflecting, boundary, but otherwise the 
results are mostly similar. In particular we conclude: 

1) Suprathermal particles can be injected very effi¬ 
ciently into the CR population via the thermal leakage 
process, so that typically a fraction of 10“"^ — 10“^ of 
the particles passed through the shock become CRs for 
0.2 < e < 0.4. 

2) For a given value of Mg, the injection efficiency 
is higher for larger values of e {i.e., weaker waves) and 
so the CR acceleration is more efficient. For example, 
in the shocks with Mg = 3, the ratio Pc, 2 /Po{u[gf be¬ 
comes 0.1, 0.16, and 0.24 for e = 0.25, 0.3, and 0.4, 
respectively. But this dependency is weaker for higher 
Mach numbers, so the acceleration efficiency becomes 
approximately independent of e in the strong shock lim¬ 
its {Mg ^ 30). 

3) For a given value of e, the acceleration efficiency 
increases with Mg, but it asymptotes to a limiting value 
for Mg ^ 30. For example, the model with e = 0.2, the 
ratio of Pc, 2 /Po{u'g)^ becomes 0.03, 0.3, 0.5, 0.8, and 
0.9, and the ratio of Pc 2 /Pg 2 approaches to 0.05, 0.6, 
1.5, 4.9, and 13 for Mg = 3, 6.8, 13.3, 40, and 133, 
respectively. 

4) In the strong shock limit of M > 30, the CR pres¬ 
sure can dominate over the gas pressure and induce a 
significant precursor where the preshock flow is decel¬ 
erated adiabatically. 

5) The CR pressure seems to approach a time asymp¬ 
totic value when a balance between acceleration/injection 
and diffusion/advection processes is achieved, resulting 
in an approximate “self-similar” flow structure. This 
is achieved in a time scale comparable to the accel¬ 
eration time scales for the mildly relativistic protons 
{p/m-pC ~ 10), which is much shorter than the cosmo¬ 
logical time scale. 

We suggest that the CR acceleration at the cosmic 
shocks are innate to collisionless shock formation pro¬ 
cess and CRs can absorb a significant fraction of dy¬ 
namical energy associated with the gravitational col¬ 
lapse during the formation of large scale structure. 
For strong accretion shocks of Mg > 10, CRs can ab¬ 


sorb most of shock kinetic energy and the accretion 
shock speed can be reduced up to 20 %, compared to 
pure gas dynamic shocks. Although the amount of ki¬ 
netic energy passed through accretion shocks is small, 
since they propagate into the low density intergalactic 
medium, they might possibly provide acceleration sites 
for ultra-high energy cosmic rays of E > 10^®eV. For 
internal merger shocks of Mg < 3 the energy transfer 
to CRs should be less than 10-20 % of the shock ki¬ 
netic energy at each shock passage, with an associated 
CR particle fraction of 10“^. Considering that ICM 
can be shocked repeatedly, however, the CRs gener¬ 
ated by merger shocks could be sufficient to explain the 
observed non-thermal signatures from ICM of galaxy 
clusters in radio, EUV and X-ray. This implies that 
the current understandings of cosmological hydrody¬ 
namic simulations could be modified by inclusion of 
this process at a quantitative level of order several tens 
of percent. 
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